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Abstract 

The 2D Euler equations is the basic example of fluid models for which a microcanical 
measure can be constructed from first principles. This measure is defined through finite- 
dimensional approximations and a limiting procedure. Creutz's algorithm is a microcanonical 
generalization of the Metropolis-Hasting algorithm (to sample Gibbs measures, in the canonical 
ensemble). We prove that Creutz's algorithm can sample finite-dimensional approximations of 
the 2D Euler microcanonical measures (incorporating fixed energy and other invariants). This 
is essential as microcanonical and canonical measures are known to be inequivalent at some 
values of energy and vorticity distribution. Creutz's algorithm is used to check predictions from 
the mean-field statistical mechanics theory of the 2D Euler equations (the Robert-Sommeria- 
Miller theory). We found full agreement with theory. Three different ways to compute the 
temperature give consistent results. Using Creutz's algorithm, a first-order phase transition 
never observed previously, and a situation of statistical ensemble inequivalence are found and 
studied. Strikingly, and contrasting usual statistical mechanics interpretations, this phase 
transition appears from a disordered phase to an ordered phase (with less symmetries) when 
energy is increased. We explain this paradox. 

1 Introduction 

Two-dimensional and geophysical flows are highly turbulent, yet embody large-scale coherent struc- 
tures such as ocean rings, jets, and large-scale vortices. Understanding how these structures appear 
and predicting their shape is a major theoretical challenge. The statistical mechanics approach 
to geophysical flows is a powerful complement to more conventional theoretical and numerical 
methods (see Ref. [T] for a recent review, or Refs. |2j [3j 0] for related approaches). In the inertial 
limit statistical equilibria describe, with only a few thermodynamical parameters, the main natural 
attractors of the dynamics. 

Recent studies in quasi-geostrophic models provide encouraging results: a model of the Great 
Red Spot of Jupiter [5], and explanations of several different phenomena: the drift properties of 
ocean rings [6j, the inertial structure of mid-basin eastward jets [6j, bistability in complex turbulent 
flows [7\, and so on. 

Generalization to more comprehensive hydrodynamical models, which can simulate gravity- 
wave dynamics and enable energy transfer through wave motion, would be interesting. Both of the 
aforementioned processes are essential in understanding the geophysical flow energy balance. How- 
ever, due to difficulties in essential theoretical parts of the statistical mechanics approach, previous 



methods describing statistical equilibria were up to now limited to the use of quasi-geostrophic, or 
simpler, models. In order to study the statistical mechanics of those models, it would be useful to 
be able to rely on numerical sampling of their microcanonical measures. 

The 2D Euler equations can be formulated in terms of the vorticity field. Points in the vorticity 
field are coupled through a long-range interaction. In contrast with traditional systems, long-range 
interaction systems are well known to show generic inequivalence between microcanonical and 
canonical ensembles ((8j El [TUl HU [12], see Ref. [8] for a classification). The microcanonical 
ensemble is richer than the canonical one, as the canonical equilibrium states form a subset of 
the microcanonical equilibrium states |12| IT3] . For these systems it is thus essential to be able to 
sample microcanonical measures instead of canonical ones. 

Creutz's algorithm |14j is a Monte-Carlo approach used to sample microcanonical measures of 
discrete spin systems |15[ 116], Whereas the Metropolis-Hasting algorithm samples Gibbs measures 
(canonical ensemble), the Creutz algorithm imposes an energy constrain and thus samples the 
microcanonical measure (microcanonical ensemble). Section [3] gives a proof of detailed balance 
and convergence to the microcanonical measure for Creutz's algorithm. Appendix A defines the 
classical concepts useful for the discussion of detailed balance. Appendix B also describes some 
improper interpretation of the algorithm that may lead to wrong results. 

The main aim of this paper is to discuss the first generalization of Creutz's algorithm to hy- 
drodynamical systems. The main novelty is the ability to deal with the statistical mechanics 
of fields rather than discrete variables. For this first work, we consider the 2D Euler equations 
and precisely define the microcanonical measures through finite-dimensional approximations and 
a limiting procedure. The generalization to the Quasi-Geostrophic model or the Vlasov equation, 
the microcanonical measure definition and their sampling through Creutz's algorithm would be 
straightforward. The method is extremely robust and could also be easily generalized to more 
complex models like the 3D axisymmetric Euler equations or the Shallow Water model. 

The 2D Euler equations and related models have an infinite number of conserved quantities, 
called Casimirs. Michel and Robert |17| have first discussed the use of large deviation theory as 
a justification of mean field variational problems for the microcanonical measure. Later on, Ellis 
and collaborators |T8] have defined approximations of equilibrium measures, through a description 
of the fields over a lattice, where A 2 is the number of degrees of freedom (lattice site). In this 
work, Ellis and Turkington have treated the energy constraint microcanonically and the Casimir 
constraints canonically. They have proven that these approximate measures verify a large deviation 
principle, where A 2 is the large deviation rate. In this paper, we study a A 2 -degrees-of-freedom 
discretized approximation of equilibrium measures (following Ellis and collaborators) , but treating 
all constraints microcanonically (as did Robert and Michel). This slightly different presentation is 
an improvement, as proceeding through discretization provides a clear and straightforward defini- 
tion of the microcanonical measures, and as the set of microcanonical equilibrium states includes 
the set of canonical equilibrium states (see the beginning of the introduction). The limit is then 
an invariant measure of the 2D Euler equations |19J . 

As was already clear in previous works |17[ [T8] (please see a detailed discussion in Ref. p]), 
the 2D Euler equations show mean-field behavior. It is therefore natural to define macrostates 
through coarse-graining of microstates. The most probable macrostate maximizes an entropy 
functional with energy constraints. The mean-field entropy for the macrosate is justified as being 
the opposite of a large-deviation rate function, where A 2 is the large-deviation rate. We explain 
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those theoretical results and their justification at a heuristic level in this paper. 

Those theoretical results (the concentration of most microstates on a single predicted macrostate 
maximizing a mean-field variational problem) provide a very interesting case for testing the Creutz's 
algorithm with numerical results. This test possibility was our main motivation for devising this 
algorithm for the 2D Euler equations, before generalizing to more complex models for which mean- 
field type large-deviation results are not yet available. 

We have checked that numerical results are in full agreement with the theoretical predictions 
from the mean-field variational problem. Independently of the mean-field variational problem, 
Creutz's algorithm provides a very simple and robust way to sample micro canonical measures. For 
instance, we have used it to describe previously unknown first-order phase transitions between 
dipole and parallel flows in a doubly periodic domain. 

Previous works considered Monte-Carlo simulations for the 2D Euler equations or related mod- 
els (see for instance a very interesting application to oceans in Ref. [20] and references therein, 
or Ref. [2]). However, those works always sampled the canonical Energy-Enstrophy measures (us- 
ing the Metropolis-Hasting algorithm and not considering other invariants). Those method thus 
provide only a very small subset of the equilibrium measures of those models. 

Dubinkina and Frank |21| recently proposed a very nice particle-mesh algorithm for the 2D 
Euler equations that conserves the vorticity distribution. This algorithm is also a way to sample 
micro canonical measures, as was shown in their paper. As a positive point, the dynamics of the 
particle-mesh method is a good approximation of the 2D Euler dynamics for finite times (whereas 
the Creutz algorithm is just a sampling of the microcanonical measure). One drawback of the 
particle-mesh algorithm is that ergodicity has to be assumed. Moreover, it seems that this particle- 
mesh approach with potential vorticity conservation has so far not been proven to be generalizable 
to more complex models (for instance the Shallow Water model). 

We also note that other numerical algorithms exist to compute equilibrium states for the 2D 
Euler equations: the Turkington and Whitaker algorithm |22| [23], relaxation equations j2^, or 
continuation algorithms 0E3J (please see Ref. [1] for a description of those algorithms). However, 
these three algorithm compute extrema or critical points of the mean-field variational problem. 
They thus rely on the large-deviation theoretical results that are so far not known to exist for more 
complex models, like for instance the 3D axisymmetric equations or the Shallow Water models. 

In this paper we also discuss the discovery of a first-order phase transition in the microcanon- 
ical ensemble, found using Creutz's algorithm for the 2D Euler equations. This phase transition 
is striking in many respects. It is a first-order phase transition in the microcanonical ensemble, 
which is a thermodynamical curiosity (see Ref. |8j). As discussed in detail in [8], such a first-order 
phase transition in the microcanonical ensemble is a sign of ensemble inequivalence, as the entropy 
curve can not be concave at such a transition point. Moreover, it is a transition from a disordered 
state towards an ordered state when energy is increased, in contrast with what could be expected 
from classical statistical mechanics arguments. This paradox is due to the negative temperature 
of the system. Indeed, then entropy, the measure of disorder, decreases when energy increases. 
We discuss this point further in section [5] From a fluid mechanics point of view, this transition 
is a drastic change of the flow topology from a dipole to a parallel flow. A very interesting recent 
example of two different phase transitions in two different statistical ensembles is discussed in Ref. 
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The outline of this paper is as follows. In Section [2] the 2D Euler equations are introduced. 
The statistical mechanics theory is treated, as well as the finite-dimensional approximation of the 
2D Euler microcanonical measure. Section [3] provides a proof of why Creutz's algorithm samples 
microcanonical measures. Theoretical predictions from the microcanonical mean-field variational 
problem presented in Section [2] are confronted with numerical results in Section |4j where we focus 
mainly on the negative temperature of the system. In Section [5] examples of phase transitions 
and an example of ensemble inequivalence are discussed. In this section, we also discuss the 
transition from a disordered state to an ordered one upon increasing energy. Section [6] provides 
some perspectives and we give comments for future work. 



2 Statistical mechanics of the 2D Euler equations 

2.1 The 2D Euler equations and invariants 

The 2D Euler equations are given by 

dtOJ + v • Vw = 0, v = e z x V*0, and oj = Aip, (1) 

where u) = (V x v) • e z is the vorticity, and v is the non-divergent velocity expressed as the curl 
of the stream function tp. The stream function is defined up to a constant, which is set to zero 
without loss of generality. The relation oj = Alp is complemented with doubly-periodic boundary 
conditions on a domain given by 3) — [0, 1) x [0, 1) and r = (x,y). The energy of the flow reads 

= \f d 2 rv 2 = \l d 2 r(VV) 2 = ~\ I d 2 r^. (2) 

The energy is conserved, i.e. dt$ = 0. The equations also conserve an infinite number of func- 
tional , named Casimirs. These are related to the degenerate structure of the infinite-dimensional 
Hamiltonian system and can be understood as invariants arising from Noether's theorem [27J. 
These functionals are of the form 

%[cu]= I d 2 rs{u), (3) 
J® 

where s is any sufficiently smooth function depending on the vorticity field. We note that on a 
doubly-periodic domain the total circulation is zero: 
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d 2 ruj = 0. (4) 



A special Casimir is 

C{a) = I d 2 rH(-u + a), (5) 



v 

where H(-) is the Heaviside step function. This Casimir returns the area of B{a) = {r | uir) < a} , 
i.e. the domain of all vorticity levels smaller or equal to a. C{a) is an invariant for any a 
and therefore any derivative of it as well. Therefore, the distribution of vorticity, defined as 
D(a) = C'(cr), where the prime denotes a derivation with respect to a, is also conserved by 
the dynamics. The expression D(a)da is the area occupied by the vorticity levels in the range 
a < oj < a + da. 
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Moreover, any Casimir can be written in the form 

C f [u] = J daf(a)D(a). (6) 

The conservation of all Casimirs (Eq. Q) is therefore equivalent to the conservation of D(a). 

The conservation of the distribution of vorticity levels can also be understood from the equations 
of motion, c.f. Eq. ([!]). We find that Doj/dt = 0, showing that the values of the vorticity field are 
Lagrangian tracers. This means that the values of u are transported through the non-divergent 
velocity field, thus keeping the distribution unchanged. The Casimirs and energy are the invariants 
of the 2D Euler equations. Their existence plays a crucial role in the dynamics of the system. 

From now on, we restrict ourselves to a if-level vorticity distribution. We make this choice for 
pedagogical reasons, but the generalization to a continuous vorticity distribution is straightforward. 
The K-level vorticity distribution is defined as 

K 

D(a) = Y,A k 6(a-a k ), (7) 

k=l 

where A k denotes the area occupied by the vorticity value a k . The areas A k are not arbitrary 
as their sum must obey Y2k=l At = \@\= 1 (the total area of the domain is equal to unity). 
Moreover, the boundary condition (Eq. Ml)) imposes the constraint ^2 k =i Ak&k = 0. 



2.2 Micro canonical measure 

As the 2D Euler equation is a conservation law (the time derivative of lv is the divergence of a cur- 
rent), it verifies a formal Liouville theorem |28( I29|. This formally justifies that the microcanonical 
measure will be dynamically invariant. In order to properly identify a microcanonical measure, we 
will thus discretize the vorticity field in finite-dimensional space, with N 2 de grees of freedom, and 
then take the limit N — > oo. As any point in the physical space has a symmetric role, a uniform 
grid has to be chosen in order to preserve the volume conservation corresponding to the Liouville 
theorem. 

We denote the lattice points by = (jj, ^J, with < i, j < N — 1 and denote co^ = 
to be the vorticity value at point r^. The total number of points is N 2 . 

As discussed in the previous section, we assume D(a) = X^fe^i J 4fc^( cr — a k)- F° r this finite- N 
approximation, our set of microstates (configuration space) is then 



X N = {uj n = (oJij)o<ij<N-i I Vi, j ujfj G {(7i, . . . , or) , and VA; # {u)$ \ = a k } = N 2 A k ) . 

(8) 

Here, ^{A) is the cardinality of set A. We note that Xn depends on D(a) through A k and a k (see 
Eq. ([7])). We note that all microstates in Xn have the proper vorticity distribution. 

In order to define a microcanonical ensemble, we also have to impose an energy constraint. 
Using the above expression we define the energy shell T^(E, AE) as 

T N (E, AE) = {u N eT N \E < g N [uJ N \ <E + AE} , (9) 

where 
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N-l N-l 
^=2^EW) 2 = -2^ E (10) 

i,J=0 i,j,i',j'=0 

is the finite- iV approximation of the system energy, with = v(r.y) being the discretized velocity 
field, AE is the width of the energy shell, and where Gij^ji is a finite— iV approximation of the 
Laplacian Green function on domain $>. We shall define these finite- N approximate fields more 
precisely in Section [3] Note that a finite width of the energy shell is necessary for our discrete 
approximation, as the cardinality of Xn is finite. Indeed, the set of accessible energies on Xn is 
also finite. Let A^E be the typical difference of between two successive achievable energies. We 
therefore assume that A^E <C AE <C Eq. 

The fundamental assumption of statistical mechanics states that all microstates in this ensemble 
are equiprobable. By virtue of this assumption, the probability to observe any microstate is 
QJ}(Eq, AE), where £In(Eq, AE) is the number of accessible microstate and is defined as the 
cardinality of the set Tn(Eq, AE). The finite- iV specific Boltzmann entropy is then given by 

S N (E ,AE) = ^lo g n N (E ,AE). (11) 

The microcanonical measure is defined through the expectation values of any observable si ' . For 
any observable s/[u] (for instance a smooth functional of the vorticity field), we refer to its finite- N 
approximation by £/n[ui n ]. The expectation value of f° r the microcanonical measure reads 

(^(^, AB), ^(/)) = (^(/)) w =_i_ Yl ^n(oj n ). (12) 

The microcanonical measure \i for the 2D Euler equation is defined as a limit of the finite-iV 
measure: 

< lx(E Q ), s/(u) >= lim < fi N (E ,AE), #J N (u N ) > . (13) 

N— >oo 

The specific Boltzmann entropy is defined as 

S(Eq) = lim S N (E ,AE). (14) 

N— >oo 

The limit measure and the entropy just defined are expected to be independent of AE in the limit 
N — > oo. This will be justified in next section with large-deviation principles. 

2.3 Sanov theorem and mean-field entropy 



Computing the Boltzmann entropy by direct evaluation of Eq. (14) is usually an intractable 



problem. However, we shall give heuristic arguments in order to show that the limit N — > oo can 



be easily evaluated. The Boltzmann entropy (Eq. (14)) can then be computed through maximizing 
a constraint variational problem (called a mean- field variational problem, see Eq. (23)). 

This variational problem is the foundation of the Robert-Sommeria-Miller (RSM) approach to 
the equilibrium statistical mechanics for the 2D Euler equations. The essential message is that 
the entropy computed from the mean-field variational problem (to be defined below) and from 



Boltzmann's entropy definition (Eq. (14)) are equal in the limit N — > oo. The ability to compute 
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the Boltzmann entropy through this type of variational problems is one of the cornerstones of 
statistical mechanics. 

Our heuristic derivation is based on the same type of combinatoric argument as the ones used 
by Boltzmann for the interpretation of his H function in the theory of relaxation to equilibrium 
of a dilute gas. This derivation doesn't use the technicalities of large-deviation theory. The aim 
is to actually obtain the large-deviation interpretation of the entropy and to provide a heuristic 
understanding using basic mathematics only. The modern mathematical proof of the relation 
between the Boltzmann entropy and the mean-field variational problem involves the theory of 
large deviations and Sanov's theorem. 

Macrostates are sets of microscopic configurations sharing similar macroscopic behavior. Our 
aim is to properly identify macrostates that fully describe the main features of the largest scales of 
2D turbulent flows and computing their probability or entropy. 

Let us first define macrostates through local coarse-graining. We divide the N x N lattice into 
(N/n) x (N/n) non-overlapping boxes each containing n 2 grid points (n is an even number, and 
iV is a multiple of n). These boxes are centered on sites = (In, Jn), where integers / and J 
verify < /, J < N/n — 1. The indices (/, J) label the boxes. 

For any microstate oj n G T^, let F^ u be the frequency to find a vorticity value er fc in box 
(I, J): 

In+n/2 Jn+n/2 

f^ n ) = ^ E E 

i=In-n/2+l j=Jn-n/2+l 

where 5d(x) is equal to one whenever x = 0, and zero otherwise. We note that for all (/, J) 

A macrostate p N = \ p¥ r T \ , is the set of all microstates oj n G Xn such 

I ' > 0<I,J<N/n-l;l<k<K 

that FN(«> N ) — PkJJ for a11 ^>^> and k ( b y abuse of notation, and for simplicity, p N — 

refers to both the set of values p? T T and to the set of microstates 
0<I,J<N/n-l;l<k<K ' J 

having the corresponding frequencies). The property FJ?jj(oj n ) = 1 imposes a local normal- 
ization constraint VI, J J2kPkU = 1- The entropy of the macrostate is defined as the logarithm 
of the number of microstates in the macrostate 

S N \p N ] = ^ log # {oj N € X N | for all /, J, and k, F kJJ (uj N ) = p^j } . (15) 

From an argument by Boltzmann (a classical exercise in statistical mechanics), using combinatorics 
and the Stirling formula, the limit JV>n>l, the asymptotical entropy of the macrostate is 

s N \p N ] A ~ i 



I— oo otherwise, 



where ^V[p^j] = J2kPklJ- ^ n large-deviation theory, this result could have been obtained using 
Sanov's theorem. 

The coarse-grained vorticity is defined as 

In+n/2 Jn+n/2 



UJjj 
1J 



72 E E <*■ (16) 

i'=In-n/2+l j'=Jn-n/2+l 
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Note that, over the macrostate p , the coarse-grained vorticity depends on p only: 

K 



UJ 



N 
13 



for 



k=l 



,iVl 



We now consider a new macrostate (p , .Eo) which is the set of microstates uj with energy (t>n[uj 
verifying Eq < $n[uj ] < Eq + AE (the intersection of Tff(E, AE) and p ). For a given macrostate 
p^, not all microstates have the same energy. Thus, the constraint on the microstate energy cannot 
be recast as a simple constraint on the macrosate p N . Therefore, treating the energy constraint 
requires a more subtle approach. The energy ( 10 ) is 



1 

'2N^ 



N-l 

E 

,j'=0 



N r* N 



(17) 



Then, in the limit JV>ii>l, the variations of Giji'j' for (i',f) running over the small box (/, J) 
are vanishingly small. Hence, Giji'ji can be well approximated by the average value over the boxes 



G 



Gijjij* + o 



n 



(18) 



From Eq. (17), using Eqs. (16 18), it is easy to conclude that in the limit iV>n>l the energy 
of any microsate of the macrosate p N is very well approximated by the energy of the coarse-grained 
vorticity 

2 N/n-1 



iVi 



n 
2N 2 



I,J=0 



where ipfj is the stream function, and is related to the velocity field and vorticity field by Eqs. (111). 



Note that in the above equation we made use of the relation 
Hence, the Boltzmann entropy of the macrostate is 



S N [p N ,E 



iV>n>l 



y N \p N ] if VA; s/ N (p%) = A k , JK\p N ] 
— oo otherwise. 



1 and 



En 



(19) 



Consider P]\f t E (p N ) to be the probability density to observe the macrostate p . By definition of 
the microcanonical ensemble and of the entropies Sn(Eq) (see Eq. (11)) and Sn(p N , Eq) (see Eq. 
(fl9l)), we have 



Pn,e (p N ) = exp {N 2 [S N [p N ,E ] - S N (E )] } . 
Let Pm be the probability density for the random variable Xm- The statement 

1 



lim 

M->oo 



M 



log [P M (X M = x)]=I(x) 



(20) 



(21) 



is called a large-deviation result. I(x) is the large-deviation rate function, and M the large- 
deviation rate. From this definition, we see that formula (20) is a large-deviation result for 



S 



macrostate p for the macrocanonical measure. The large-deviation rate is N and the large- 

?o). 



deviation rate function is —Sjy\p , Eq] + Sn(Eq 



We now consider the continuous limit n — > oo, N — > oo. The macrostates p^ are now seen as 
finite- N approximation of p^, the local probability to observe w(r) = : Pfc(r) = (5(u)(r) — a^))- 
The macrostate is now characterized by p = {pi, . . . ,Pk}- Taking the limit JV>ti> 1 allows us 
to define the entropy of the macrostate (p, Eq) as 

cr r,i_l y M= -EkUdrPklogPkif Vk^\p k } = l, ^ (p k ) = A k and gp]=E 
I — oo otherwise, 

where Vr, ,jV(y) = ^2 k= iPk{ r ) = 1 is the local normalization. In the same limit, it is clearly seen 



from definition (15) and result (22) that there is a concentration of microstates close to the most 
probable macrostate: the equilibrium state. The exponential concentration close to the equilibrium 
state is a large-deviation result, where the entropy appears as the opposite of a large-deviation rate 
function (up to a constant). 

The exponential convergence towards this most probable state also justifies the approximation 



of the above entropy with the entropy of the most probable macrostate, Eq. (14), as 

S(E ) = max f S*\p] | S[u)} = Eq, Vfc £/{p k ) = A k \ , (23) 

{p}\._,V (r)=l I J 

where p = {p\, . . . ,Pk}, Vr ,jV(y) = ^A^i^( r ) = 1 i s the local normalization, y[p] is as defined 



in Eq. (22), and £/[pk] is the area of the domain corresponding to the vorticity value U) = a^. The 



fact that the Boltzmann entropy S(Eq) (Eq. (14)) can be computed from the variational problem 



( 23 ) is a powerful non-trivial result from large-deviation theory. 

In the next section, we shall test the prediction of concentration of microstates close to the 
equilibrium macrostate with numerical simulations. We first define Creutz's algorithm and explain 
why it is able to sample microcanonical measures. We continue by applying this algorithm to the 
2D Euler equations. 



3 Creutz's algorithm 

Creutz's algorithm was introduced by M. Creutz in 1983 [14j. It is a generalization of the 
Metropolis-Hastings algorithm that samples the Gibbs measure (with a Boltzmann factor). Creutz's 
algorithm, on the other hand, samples the microcanonical measure in the energy shell Eq < <§ < 
Eq + AE (a uniform distribution over the set Tn(Eq, AE)). Here, the system energy is denoted 
by S. 



In Section 3.1 we present Creutz's algorithm and prove that it actually samples the micro- 



canonical measure. In Section 3.2 we provide a method to calculate the inverse temperature f3 
using this algorithm. We shall refer to Appendix A for a precise definition of Markov chains, the 
detailed balance condition, and invariant measures (stationary distributions). 

Creutz and others using his algorithm used the notion of a daemon. The aim of the 

daemon is to allow for slight energy fluctuations, which are necessary for systems with a discrete 
configuration space. With our notation the daemon energy E& is nothing else than E<± = Eq + 
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AE — £. The original Creutz algorithm samples a uniform measure over all microstates of energy 
smaller than Eq + AE. For this measure, in systems with positive temperature and a large number 
of degrees of freedom, the energy distribution is concentrated close to Eq + AE, and typical energy 
fluctuations are small. 

In the 2D Euler case considered in this paper the temperature can be negative, as will be 
discuss in Section 4,1,2| Microstates with the smallest possible energy then become overwhelmingly 
probable. In order to sample the microcanonical measure, we then need to impose a lower-bound 
on the energy. In order to cope with all possible temperature cases we sample a uniform measure 
over the energy shell Eq < £ < Eq + AE. Because of energy concentration properties, the 
microcanonical measure is independent on those definitions or on the value of AE in the limit of 
an infinite number of degrees of freedom. 

Classic heuristic arguments using the daemon energy lead to misleading conclusions, for instance 
in the case of negative-temperature systems. For this reason we prefer not to use this concept at 
all. We notice moreover that the deamon concept is actually not useful as its energy contains no 
more information than the state energy £,. 



3.1 Definition of Creutz's algorithm 

We consider an ensemble of states of a physical system. Each state x is a set of M values: 
x = (xj)i<j<M- The set of states 



X = {x = (x f )i<<<M} (24) 

is called the configuration space. Note that i £ N is the index for components of x. 

For instance, a state of an Ising spin system is given by xi £ { — |, + 2}- Note that the values 
x could also be continuous, i.e. \/i Xi £ R. For the 2D Euler equations, the configuration space is 
X N as defined in Eq. (|]), and M = N 2 . 

The energy of each microstate is given by $(x). Furthermore, we define 



T(E , AE) = {x £ X I Eq < S{x) < E + AE} (25) 

as the set of microstates in the energy shell Eq < <§{x) < Eq + AE. The microcanonical measure 
is defined as the uniform measure over V(Eq, AE) (each microstate in T(Eq, AE) has probability 
Q(E, AE) = (#T(Eq, AE))^ 1 to occur, where we recall that if=(A) is the cardinality of set A). 
Furthermore, as in the previous section, we assume that AmE <C AE <C Eq. Our goal is to sample 
the microcanonical measure over T(Eq, AE). 

We assume that there exists a Markov chain defined through a sequence of random numbers 
{y l £ X} 1>() with transition probability T(x,x') (so T(x,x') is the probability to observe y l+1 = x 

if y l = x'). Here, I £ N is the index for the position in the Markov chain. We assume that 
3~ verifies detailed balance for the uniform measure on X. This is equivalent to the statement 
T(x,x') = T{x\x). See Appendix A for more details. Lastly, we note that 3F (or T) does not 
depend on the system energy. 

In order to sample the microcanonical measure over T(Eq, AE) an algorithm is needed that gen- 
erates realizations { z l } of a new Markov chain J3 defined by its corresponding transition probability 
Q and over the configuration space T(Eq, AE). We shall call this algorithm Creutz's algorithm, 
and we define it in the following way. 
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Let the system's current state be y l = x'. We pick at random x € X with probability T(x,x'). 
If Eq < ${x) < Eq + AE then we accept this move and y l+1 = x. Otherwise, we do not accept 
and y l+l = y l = x'. Iteration of this procedure defines a Markov chain J2. 

We can show that Q verifies detailed balance for the uniform measure over T(Eq, AE). It is 
easily checked that Q(x, x') = T(x, x') if x ^ x', and Q(x, x) = T(x, x) + J2 x 'ex\r(E AE) , %) = 
1 — J2 x 'er(E AE) xj^x' x ')- (Q( x i x ) i s the probability of a move from y l to y l+1 to fail because 
of the energy constraint.) We thus find that \/x,x' £ T(Eq, AE) : Q(x,x') = Q(x',x) by virtue 
of the detailed balance on 3? . Therefore, verifies a detailed balance for a uniform measure 
on T(Eq, AE). Thus, has this uniform measure as a stationary measure, see appendix A. In 
conclusion, if £} is ergodic then it samples the microcanonical measure over T(Eq, AE). 

We conclude by describing how to properly empirically sample an observable A. The expecta- 
tion value of observable A is computed through 



<A(y)>= lim jJ2 A (y 1 )- ( 26 ) 

1=1 

It is important to notice that if for y l , the Creutz algorithm fails to change the state K times 
(y = U l+1 = ■ ■ ■ y l+K ~ 1 ), then these K configurations need to be included into the sum (26). We 
clarify this statement in Appendix B. 

3.2 Temperature computation using the Creutz algorithm 

We continue with the computation of the inverse temperature from the energy distribution at 
equilibrium. We denote pm(E) as the density of states for energy E. The number of microstates 
in T(Eq, AE) then reads 

pEo+AE 

n M (E ,AE)= / dE pm(E). (27) 

JEo 

We assume that the energy density of state has a large deviation behavior, i.e., 

VE lim ~logp M (E) -> S(E). (28) 

M->oo IVl M^foa 



We can assume a stronger property, namely 

,MS{E) 

We note that this assumption is verified for most microcanonical measures. For the 2D Euler 



p M (E) ~ C(E)e M ^ E >. (29) 

M— >oo 



equations, it follows from the large-deviation results discussed in Section 2.3 

In the microcanonical ensemble the temperature is defined as /3 = §^ It then follows that 
S(E) = S{E ) + P(E - E ) + o(E - E Q ) and p M (E) ~ C ( E } e MS(E )+pM(E-E ) We thug hav6) 

M— >oo 



using the fact that Eq. (27) is a Laplace integral, that Huim^-oo log £Im(Eq, AE) = S(Eq) is 
independent of AS. 
Moreover, 

p M (E) ~ C 1 (E ,AE)e-P M ( E - E °\ (30) 
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and 



P(E) ~ C 2 (E ,AE)e-P M ( E - Eo l (31) 

M— >oo 

The energy distribution is exponential with rate (3M. If /3 is positive, the energy of the system is 
concentrated close to Eq + AE. If /3 is negative, the energy of the system is concentrated close to 
Eq. 

We will return to this result in Section |U 
3.3 2D Euler algorithm 

We now want to apply the Creutz algorithm to the 2D Euler model. We will follow the general 



definition from Section 3.1 Thus, we need to (i) define a Markov chain for the 2D Euler model 
and (ii) precisely define the approximate energy such that we can sample states in the energy shell 
[E ,E + AE]. 

Firstly, following the previous section, we consider the Markov chain 3T, now defined through 
its configuration space 



X N = {u N = (uij)o<i,j<N-l I Vi, j ufj £ {<Ti, ... , a n } , andV/c # {ujfj \ ujf- = a k ] = N 2 A k ] , 

(32) 

and transition probability T(oj,oj'). T denotes the probability to go from a vorticity state y l = 
uj' £ Xn to y l+1 = uj £ Xn, where as before I > denotes the position of y l in . We assume the 
detailed balance condition holds for ST . 



Secondly, we need to calculate the approximate energy $n given by Eq. (10). The discretized 
vorticity field uj n is transformed to Fourier space with a two-dimensional (discrete) Fast Fourier 
Transform (FFT). The discretized velocity field and stream function field are then computed by 
making use of the Fourier representation of Eqs. ([!]) and an inverse FFT, resulting in v N = 

( v ij)o<iJ<N-i and ip N = (tp(j)o<i,j<N-i- 

With these definitions the microcanonical measure over Eq. (§ can be sampled. We take the 
current system state to be y l = uj £ Xn (the current vorticity field configuration) and y l+1 = uj £ 
X/v as the next state (next vorticity field configuration) in ST. Similarly, the current energy is 
denoted by <£jv[z/ ] an d the energy of the next state by <£at[?/ +1 ]. 

We now describe how we construct a realization of the Markov chain under the energy 
constraint. The vorticity field y° = uj n £ Xiq is initialized with N 2 /2 values of cr + = +1 and 
N 2 /2 values of <r_ = —1 such that conditions Ylk=i = 1 an< i Sfj=o ^ij = ® are satisfied (see 



Section 2.1). We use a two-level potential vorticity distribution as an example, but generalization 
to higher-level distributions is straightforward. To sample the microcanonical measure, obeying 
the energy and vorticity distribution constraints, the following routine is followed: 

1. Let the current state of the system be y l = uj' . We randomly chose two lattice sites, (io,jo) an d 

ji) with a uniform distribution. The values at these sites are u'i j and w^^, respectively. 
The vorticity values at these positions are then interchanged. In other words: ujij = uj'- 
Vij / Oojo, nil}, ^ioj = w- Ul , and uj hjl = uj' ioj(j . We note that by interchanging two 
vorticity values, the new field uj still belongs to X^. 

2. The new energy, (on[oj n ], is calculated. 
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3. The energy check is performed. If Eq < S^[uj n \ < Eq + AE then we accept the move in 
step 1, i.e., u! l+1 = u. If the new energy is not in the allowed energy range, we reverse step 
1 such that we get the old vorticity field back, i.e., = ur = uj\ The new energy in this 
case is thus given by <^v"M = $n{u'\- This step ensures conservation of energy and only 
configurations are allowed which are in the set Fn(Eq, AE), see Eq. d9J). In either case, we 
return to step 1. 

Iteration of steps (1-3) then builds a realization {u/}^ >0 of the Markov chain Q. Since we assumed 



that detailed balance holds for 2? ' , we know (by virtue of Section 3.1 ) that the microcanonical mea- 
sure is an invariant measure of Q. If we assume Q to be ergodic then we sample the microcanonical 
measure. 

This vorticity exchange method has some analogies with spin ones in the Kawasaki algorithm 
[30] . The Kawasaki algorithm is a convenient way to treat dynamics with conservation laws in 
spin or lattice gas Monte Carlo dynamics. Our algorithm slightly differs from Kawasaki's, as we 
consider non-local vorticity value exchanges. 

We define a Monte-Carlo step as N 2 accepted changes in the vorticity field, for reasons men- 
tioned in Section [3.1 1 Typically, equilibrium is reached after ten Monte-Carlo steps for a two-level 
distribution and fifty for a three-level vorticity distribution. We note that the equilibration time 
depends on the energy Eq chosen. 



4 Numerical results 

In this section we present the numerical results obtained by sampling the microcanonical measure 
of the 2D Euler equations using the algorithm described in the previous section. For pedagogical 
reasons, we restrict the simulations to two- and three- level vorticity distributions, but note that 
generalization to more general vorticity level distributions is straightforward. We consider the 
vorticity distribution (see Eq. ([7])) 

D{a)= l -5{a-a + )+ l -5{a-a-) (33) 
for the two-level distribution, and 

D(a) = \5{a - a+) + \s(a - a ) + \d{a - a-) (34) 
for the three-level distribution. In the above, a + = 1, cr = 0, cr_ = —1. 



4.1 Mean-field behavior and negative temperature 
4.1.1 Mean- field predictions 



The variational problems (23) associated to the two- and three-level cases have been explicitly 
solved in the past, see Refs. |24[ [25j [31]. We summarize their results here and compare them to 
our numerical results. 

For the two-level vorticity case the local probability to find a vorticity value of to = a+ = 1 is 
found to be p\ = 1(1 — tanh(a — /3ip(r))) [5], where a and /3 are Lagrange parameters associated 
with the conservation of area A and energy, respectively. Furthermore, f3 is the inverse temperature: 
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P = Since there are only two possible values for oj, the local probability to find a vorticity 
value of oj = er_ = — 1 is therefore p2 = 1 — p%. Hence, the locally averaged vorticity reads 

oj = p\a + + (1 — = — tanh(a — Pip). (35) 

We consider a symmetric distribution D(a) (D(a) = —D {—&)). If we assume this symmetry then 
a = 0. 

For the three-level distribution a similar calculation, assuming symmetry is not broken, yields 

m 

_ _ /xsinh(^) 
UJ ~ l + M cosh(/3V0' 
where /x is an additional Lagrange parameters arising from the conservation of A. 

4.1.2 Negative temperature in doubly periodic domains 

The first statistical mechanics approach to self organization of 2D turbulence was Onsager's work 
on the point vortex model |32j . Onsager argued in his paper that ensembles of point vortices may 
have a negative temperature. This property is traced back to the fact that, in contrast with many 
other systems, the phase space is bounded (the point vortex model has no quantity analogous to 
kinetic energy in particle systems, allowing the system to explore higher and higher energy with 
increasing entropy). The same argument also holds for the 2D Euler equations with continuous 
vorticity fields, so one expects negative temperatures to exist for some ranges of energy. In this 
section we show that, in the case of a doubly-periodic domain, the inverse temperature of the 2D 
Euler equations is actually always negative. Note that the proof has first been established by A. 
Mikelic [33] • 

From the mean-field variational problem we can prove that for any vorticity distribution there 
exist functions / such that 

u = A$ = f(P1>). (37) 

Moreover, it can be proven that f'(x) > (see for instance Ref. [5]). 

Multiplying the above equation by Aip and integrating by parts gives 

Using f'{x) > 0, we conclude that P < 

Let us find the value of P in the low energy limit for domain @. Previous works [lj have 
shown that in this limit, the expression for the vorticity field for the parallel flow is given by 
oj = Acos(2irx + (p) (or oj = Acos{2ny + (f))), where A is a constant depending on E and <p is 



an arbitrary phase. Furthermore, in the same limit, we may linearize Eq. (37) and find oj ~ Pip. 
From oj = t4cos(27tx + </>) = Aip we compute ip = —(Att' 2 )~ 1 oj. Together with oj ~ Pip we find that 
P = —Air 2 , which is indeed negative for all energies. 



1 Note that doubly-periodicity is required to ensure the temperature is always negative. This is due to a partial 
integration step in the proof. 
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4.1.3 Numerical results 



We now confront the analytical results with numerical computations. In order to do so, we must 
first ensure that the system is in equilibrium. We can test this by computing the mean-field 
entropy. We recall the entropy of a coarse-grained macrostate p N , see Eq. (19). The mean-field 
entropy is computed for the two-level distribution [K = 2,a+ = l,cr_ = —1) with parameter 
values N = 256 and n = 3, and energies ranging from Eq = 0.01E max to Eq = 0.99E max . Here, 
Emax is the maximum possible value of energy the two-level system can obtain. We will compute 
Emax below. We take AE = 0.01-Eb but note that its value is not of high importance since the 
actual energy fluctuations at this resolution are much lower than AE. 

The result is depicted in Fig. [ija). Equilibrium is reached quickly for all energy levels. To 
be on the safe side, we take an equilibration time of one-hundred Monte-Carlo steps. We define 
a Monte-Carlo step as M accepted iterations of the Creutz algorithm. With this definition, the 
size of a state x and the Monte-Carlo step scale linearly, such that after K Monte-Carlo steps each 
value Xi has been changed K times on average. All results that follow below have been obtained 
after the system has reached equilibrium. 
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Figure 1: (a) Mean-field entropy as a function of Monte-Carlo steps for a two-level vorticity 
distribution at a resolution of N = 256. Equilibrium is reached within fifty Monte-Carlo steps for 
any energy level, (b) Mean-field entropy as a function of Eo/E maa: . The numerical fluctuations are 
smaller than the size of the balls. 
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Two-level vorticity distribution. For the 2D Euler equations and for a given vorticity distri- 
bution, energy maxima correspond to segregated states. For example, when the left half of the 
vorticity field contains only values of <r+ = 1 and the right half consists solely values of <r_i = — 1 
we have Eq = E max . 

We calculate E max analytically. We have cPtp/dx 2 = — 1 for < x < \ and d 2 tp/d 2 x = +1 for 
i < x < 1. The general solutions to these two differential equations on their respective domains 
are denoted by and ip\. These equations are complemented with the boundary conditions 



(x) 



x 
' 2 



+ tI and 



l 



The numerical 



V^-i(O) = Vi(l) = 0, = ^l(i), and = ^(i). We find ^_ 

ipxix) = \ — \x + -r. The maximum energy is now easily calculated: E mc 
value of E max = 0.0104 corresponds well to this theoretical value. 

We perform numerical simulations at N = 256 for Eq = 0.9Emo,x an d an energy tolerance set to 
AE = 0.01-Eo- After reaching equilibrium, we compute the p N and stream function i() N by point- 
wise averaging one-hundred fields which are separated in time by iV 2 permutations, respectively. 
The averaged coarse-grained vorticity field H) N is then computed from p N . Fig. [2|a) shows the 
averaged coarse-grained vorticity Q N versus the averaged stream function ip N . 

The data points (blue) are fitted with Eq. (35) by tuning the Lagrange parameter /3. The 
result (red curve) corresponds to a temperature of (3 = —96.2. The numerical results are clearly 
in agreement with the theoretical predictions. In Fig. [2^b) the stream function field is plotted, 
showing a unidirectional flow. 
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Figure 2: Numerical result of a two-level vorticity distribution with resolution 256 x 256 at energy 
Eq = 0.90£' max and coarse-graining parameter n = 3. (a) Averaged coarse-grained vorticity versus 
stream function. An inverse temperature of (3 = —96.2 is found, (b) Averaged stream function 
field shows a parallel (pure) flow. 



The fluctuations in uj visible in Fig. are a result of the coarse graining of u (or equiva- 
lently, in p). The level of fluctuations in p is of o(l/Vn?) = o(l/n). This amounts to fluctuations 
of about 0.3 for n = 3. This is close to the observed fluctuations near the center (■*/> = 0), see Fig. 




Three-level vorticity distribution. Although determining E max for the two-level case was pos- 
sible due to the simple geometry, this is less trivial for the three-level case. Instead of trying to 
compute E max analytically, we determine it numerically. We let the algorithm run for some time 
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and accept only moves which increase the energy. We find that the system converges towards a 
maximum energy of E max = 0.005208. 

We now show results for a three-level distribution with N = 256, Eq = 0.9E max and an energy 
tolerance of AE = O.Oli^o- We follow the same averaging procedure as in the two-level case, see 
above. 

In Fig. Ma) the averaged coarse-grained vorticity co N is plotted against the averaged stream 
function tp. The data points (blue) are fitted with Eq. (36), see Fig. [3] for more details. This is 
again in agreement with theory. In Fig. [3^b) the averaged stream function field is plotted, showing 
a dipole flow. 




Figure 3: Numerical result of a three-level vorticity distribution with resolution 256 x 256 at energy 
Eq = 0.9E max and coarse-graining parameter n = 3. (a) Averaged coarse-grained vorticity versus 
stream function. Fitting parameters are f3 = —252.1 and \x = 0.0969. (b) The averaged stream 
function field shows a dipole flow. 



4.2 Inverse temperature computation 



The 2D Euler equations have the property of negative temperature as discussed in Section 4.1.2 



In this section we present three different methods to compute the inverse temperature of the 
two-level vorticity system. The first method is a direct measure of energy fluctuations, as discussed 
in Section 4.2.1 It is independent of any assumption (mean-field, large N, etc). It is therefore 
a good method to verify the two other methods which rely on mean-field approximations. The 



second method, studied in Section 4.2.2, uses the mean-field equation (35) and a fit to compute 



/3. The third method, proposed in Section [4.2.3 uses the mean-field approximation of the entropy 
([22]) to find (3 from the relation j3 = dS/dE. 

We show that all three methods give similar results. Furthermore, we find negative temperatures 



for all energies, in agreement with our finding in Section 4.1.2 



Lastly, in the high energy limit, we find that the inverse temperature tends to j3 = —00. In 
the low energy limit, we find j3 = —39.74, which is close to the theoretically predicted value of 
(3 = —Air 2 for our domain 0, see previous section. 

4.2.1 Temperature from the energy distribution in the Creutz algorithm 



This method makes use of the theory discussed in Section 3.2 and is applied to the two-level case. 
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After equilibrium is reached, the system's energy is computed ten times per Monte-Carlo step for 
a total of one-hundred Monte-Carlo steps. The histogram is shown in Fig. Qa) for Eq = 0.95E max . 
The data is fitted with Eq. (30), from which we obtain a value of f3. 

This method is repeated for different values of energy Eq, such that we obtain a plot of (3 versus 
energy, shown in Fig. (4^b). Notice that the temperature is negative over the full range of energy 
values, as expected. 



/J(£ /£ max ) 




0.009850 0.009851 O.O09B52 0009833 0.0098S4 0.009855 0.009856 



(a) (b) 

Figure 4: Two-level vorticity distribution with N = 256. (a) Energy distribution after reaching 
equilibrium. The energy distribution is proportional to Q,n(E). The distribution follows an ex- 
ponential law (Eq. (|30[)), see red curve. This gives access to the inverse temperature (3. (b) The 
inverse temperature, obtained using the method of (a), is plotted for several energies Eq. The error 
bars show the standard error of the estimate of (3 obtained via the fitting procedure. 



4.2.2 Temperature from the mean-field equation 



We compute the inverse temperature through the mean-field equation (16) by using the same 
method discussed in Section 4.1.1 see also Fig. [2ja). 

After the system has reached equilibrium, the average fields ip N and Co N are computed for 
several energies Eq. We plot Gj n versus tf) N for each energy and fit the result with Eq. (35) from 
which we obtain a negative temperature /3(E), see the red curve in Fig. [5] 
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Figure 5: Inverse temperature plotted against energy from the equilibrium energy distribution 
(blue), the mean-field equation (red) and the mean-field entropy (black). 



4.2.3 Temperature from the mean-field entropy 



This method relies on the mean-field entropy, see Eq. (22) and Section 4.1.1 

For each chosen energy Eq the algorithm is run until equilibrium is reached. The simulation 
is then continued for another fifty Monte-Carlo steps. We average the coarse-grained entropy and 
energy over this interval and repeat this process for energies ranging from Eq = to Eq = E max . 
The data is fitted with a polynomial, which can be derived. The temperature of the system is then 
found via /3(E) = ^(E). The resulting values of (3 for are then plotted against energy, see black 
curve in Fig. [5] 

From the same figure we can conclude that the temperature is indeed negative for all values of 
Eq. Furthermore, note that all three methods give very similar results, especially in the low energy 
limit. In this limit, statistical mechanics theory predicts a temperature of (3 = —Air 2 , see Section 



4.1.2 The values in this limit, shown in the figure, indeed approach this theoretical value. 



5 Phase transitions and statistical ensemble inequivalence 

In this section, we use Creutz's algorithm in order to study the microcanonical measures for the 
2D Euler equations. We are specifically interested in the study of phase transitions, as they play 
a major role in the dynamics of equilibrium and non-equilibrium flows [I]. The numerical 
computations shown in this section are compared with the low-energy statistical equilibria and 
phase diagrams analytically computed in Ref. 0H]. Creutz's algorithm allows to search beyond 
the low energy limit. In Section |5.2| we observe a first-order phase transition that does not exist 
in the low-energy limit, for a three-level vorticity distribution. 



5.1 Summary of theoretical results 

We begin with a summary of the theoretical results describing phase transitions for the 2D Euler 
equations in the low energy limit (TJ [33]. In a domain with doubly periodic boundary conditions 
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and in the low energy limit, statistical equilibria are well approximated by largest scale eigenmodes 
of the Laplacian 



u(x,y) ~ A cos(27rx + 0) + B cos(2vry + </>'), (39) 

where A, B are constants. 

For a square geometry, there are three possible equilibrium flows. Two of these equilibria have 
amplitudes of (A = 0, B ^ 0) and (A ^ 0, B = 0), respectively, and are called pure states, or 
parallel flows. The last type of flow, called a symmetric dipole, is the case for which A = B. 
Symmetric dipoles and of parallel flows have been found numerically with our algorithm, see 
previous sections. We proceed now to a more detailed study of the phase transitions between those 
equilibria when energy is changed. 

It is shown in Refs. [TJ |33] that the selection of either a dipole or a parallel flow in the low- 
energy limit is related to inflection points of the relation between the coarse-grained vorticity and 
the stream function. To be more precise, consider the Taylor expansion 

co = f(l3ifj) with f(x) = x + a 4 x 3 + o(x 3 ). (40) 



Using (3 < (see Section 4.1.2), we note that when 0,4 > 0, the curve uj — ip bends upwards for 
positive x, similar to a hyperbolic sine. When 04 < 0, then it bends downward for positive x, 
similar to a hyperbolic tangent. The theory predicts parallel flows when 04 < 0, and dipole flows 
when 04 > 0. 

Let us study this criteria in the case of a two-level vorticity distribution. From the theoretical 



predictions (Eq. (35)) we find for the two- level case with symmetric vorticity distribution (a = 0) 
that uj — Pip — |/3^i/> 3 , from which we find that 04 = —| < 0, in accordance with the tanh-like 
behavior observed (see Fig. [2^a)). We thus expect to always observe parallel flows for the two-level 
case, and no phase transitions. This is in agreement with the results obtained by using Creutz's 
algorithm. 



The three- level case is more interesting. Linearization of Eq. (36) yields uj ~ j^fitfj + 



^(i+m)^ 3 ^ 3 ' ^ or ^ < 1/2 we find that 04 > and expect to observe a dipole. This is the 
case studied in previous section (see Figs. 3). Interestingly, in the three-level case, a^/J,) = gWqz^pf 
is negative for either \i > 1/2 or // < 0. This open the possibility for a phase transition in the 
three-level case. 

We have not tried to theoretically compute /jasa function of the energy as this would be very 
involved, but we find that a phase transition actually exists using Creutz's algorithm. We show 
the result in next section. 



The theoretical results [1, 34j are valid in the low-energy limit. In this limit the theoretical 
results indicate that the transition is of second-order (symmetry breaking, by transition from 
parallel to dipole flows) and occurs when 04 = 0. In the following we perform computations both 
in the low-energy limit and for large energies. We find that the sign of 04 remains relevant for 
finding phase transitions. However, we will see that for large energies the transition can be of 
first-order. 
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5.2 Phase transition in the three-level case 



In order to study phase transition between dipole and parallel flows we first define an appro- 
priate order parameter. We consider the Fourier coefficients of the vorticity field U) , denoted 



Wk, where k = (k x ,k y ), such that in Eq. (39) A = \oo(2n,o)\ an d B = |w(o,27r)|- We define 
LOmin = min{|d>( 2vri o)| , |w( ,27t)|} and Umax = max { | a>( 27r ,o) | ! |w(o,2tt)|} ■ Tnen tne order parame- 
ter 

r = ^ (41) 

Umax 

characterizes the equilibrium state in which the system resides. We have < r < 1, and for r = 1, 

u m in = u max , A = B and we observe a purely symmetric dipole, whereas for r = 0, uj m i n = 0, we 

find that either A or B is equal to zero, corresponding to either a horizontal or vertical parallel 

flow. For an intermediate value of r, a mixed state between a parallel flow and dipole will be 

observed (which is not a statistical equilibrium but can be observed transiently). 

For N = 128 we compute the order parameter 256 times per Monte-Carlo step for energies 

ranging from Eq = to Eq = E max for the two- and three- level cases, see Fig. [4] For both 

cases we found degeneracy at very low energies. Indeed, in a square box, theory [TJ [1] predicts a 

second-order phase transition at E = as — > 0. Then for finite N, both parallel and dipole 

E— >0 

flows are possible at non-zero but small energies due to finite-size effects. Fig |4^a) confirms this. 
It is also expected in Ref. [7J [T] that there is no other phase transition for E > in the case of 
two levels of vorticity. Our numerical results agree with this observation. 

For the three-level case, Fig. p[b) shows that a phase transition exists around Eq = 0.7E max . 



In order to determine the order of the phase transition, we ran Creutz's algorithm with three 
levels of vorticity, increasing energy adiabatically from E m = 0.60E max up to Em = 0-82E max ; see 
Fig. [7] We took a resolution of N = 256 such that fluctuations are small. 




(a) (b) 

Figure 6: Order parameter r versus E/E max = to E/E max = 1.0 for a two-level (a) and three- 
level (b) vorticity distribution. Both distributions show degeneracy at low energies. Only the 
three- level case exhibits a phase transition around E = 0.7E max . 



Two simulations were run at this resolution; a forward (low to high energy) and a backward 
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(high to low energy) simulations. The energy is adiabatically increased (decreased) with a rate 
corresponding to 4.4 x 10~ 4 E max per Monte-Carlo step, for a total of 500 Monte-Carlo steps. 
Within each Monte-Carlo step, we compute averages of the order parameter over 512 realizations, 
for each of the 500 energy levels. The result is shown in Fig. [7] We observe hysteresis behavior, a 
typical signature of first-order phase transitions. 



r 




max 



Figure 7: Hysteresis in the transition from parallel flow to dipole. The plot shows the order 
parameter versus energy for a forward (low to high energy) simulation (blue) and a backward 
(high to low energy) simulation (red) (resolution N=256). 



Using the standard terminology in describing phase transitions, the dipole is the ordered phase 
in the sense that it has lost the translational symmetry. The parallel flow is the corresponding 
disordered phase. An interesting remark is that we observe a transition from the disordered phase 
to the ordered one as the energy is increased. This is in contrast with classical thermodynamics 
and statistical physics results. For instance, the first-order solid-liquid transition is a transition 
from the ordered phase (solid with broken symmetry) to disordered phase (liquid) when energy (or 
temperature) is increased. Similarly, ferromagnetic-paramagnetic phase transition is a second-order 
phase transition from the ordered phase (ferromagnetism, with broken symmetry) to a disordered 
phase (paramagnetism) when energy is increased. 

This paradox is due to the fact that our system has a negative temperature. Then entropy 
decreases with energy. It is thus natural to expect a transition from the disordered state (high 
entropy and symmetric) to the disordered state (low entropy and with broken symmetry) when 
energy is increased. This paradoxical property is traced back to the fact that, in contrast with any 
other systems, for the 2D Euler equations the phase space is bounded (there is no equivalent of 
kinetic energy allowing the system to explore higher and higher energy with increasing entropy). 

We have described a first-order phase transition in the microcanonical ensemble with hys- 
teretical behavior. In systems with short-range interactions, microcanonical first-order transitions 
usually do not exist because of the possibility of phase coexistence [8j, whereas it is a generic feature 
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in systems like the 2D Euler equations which has long-range interactions. A very general argu- 
ment [8] show that such a first-order phase transition in the microcanonical ensemble is necessarily 
associated with a situation of inequivalence between the microcanonical and canonical ensemble. 

6 Conclusions 

In this paper we have presented a novel numerical method based on Creutz's algorithm to sample 
microcanonical measures of hydro dynamical systems. Although we have only presented numerical 
results for the 2D Euler model, we stress that this numerical scheme can easily be generalized to 
more complex hydrodynamical models such as the axisymmetric 3D Euler equations or the Shallow 
Water equations. 

For the 2D Euler model, we have reproduced the (theoretically) well-known equilibrium states 
characterized by parallel and dipole flows. Using our algorithm, we were able to compute the tem- 
perature of the system in three different ways. One of these approaches allowed the computation 
of the temperature without making use of any mean-field assumptions, and therefore served as a 
verification tool for the mean-field approximations made in theoretical predictions. All three meth- 
ods match very well, showing consistency between the numerical approach and the mathematical 
results, and also proves that the mean-field description is exact in the large- N limit. 

Furthermore, we have found a previously unknown phase transition in the microcanonical 
description of the 2D Euler equations. We have shown that in the energy range where the transition 
occurs there is an ensemble inequivalence between the microcanonical and canonical ensemble. This 
transition is very interesting from a statistical mechanics point of view, as it is a transition from 
an ordered phase with broken symmetry towards an asymmetric ordered phase when energy is 
increased. From a fluid mechanics point of view, it is also very interesting, as it describes a 
discontinuous change of the flow topology. A similar phase transition was already observed in a 
non-equilibrium framework for the 2D stochastic Navies-Stokes equations [7j. This new equilibrium 
result will probably very useful in explaining why the 2D-Navier-Stokes non-equilibrium phase 
transition has the phenomenology of a first-order phase transition rather than of a second-order 
one. 

The numerical method we propose is extremely easily implemented and the generalization 
to similar models is rather straightforward. We guess that it will have many applications for 
theoretical, experimental and geophysical flows. 

Appendix A Markov chains, detailed balance and invariant distri- 
butions 

In this appendix we formalize the notions of Markov chain, detailed balance, and reversible Markov 
chain. These definitions are used in Section 13.11 and 13.31 

Markov chain A Markov chain is a (mathematical or physical) system that undergoes changes 
from one state to another between a finite number of states. Such changes are called transitions and 
the probabilities associated with the state changes are called transition probabilities. The system 
is driven by a random and memoryless process, i.e., the next state in the chain depends only on 
the current state and not on the sequence of states that preceded it. The memoryless property of 
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the system is usually referred to as the Markov property. The set of all possible states, called the 
configuration space, and transition probabilities fully characterizes a Markov chain. 
The configuration space is defined by 



X = {x= {xi)i<i<M} . 

A state is thus described by a set of M values: x = (xi)i<i<M- 

Formally, a Markov chain & is a sequence of random variables {y l € X}^ >Q with the Markov 
property, where the transition between states are determined by transition probabilities 

T(x l+1 , x\ x l ~\ ...,x°) = T(x l+ \x l ) = T(x, x). 

It gives the probability to go from a state y l = x' to the state y l+l = x. We call a sequence {z l } l>0 
a realization of the Markov chain & . Note that i G N is the index for components of x, while 
I G N is the index for the position in the Markov chain. Furthermore, we denote P l {x') as the 
probability for y l = x' . The probability to observe the system in state y l+l = x is then given by 
P l+1 (x) = Ylx'eX T(x, x')P\x'). The stationary distribution, denoted is defined such that 
P°°{x) = , }2 x i ( zx r r{ x -, x ')P' X3 { x ')- The stationary distribution is thus invariant under T. 

Detailed balance A Markov chain ST is said to be reversible with respect to the distribution 
P°°(x) if 

Vx, x' G X : P°°(x)T{x, x') = P 00 (x , )T(x', x). (42) 



This condition is also known as the detailed balance condition. Summing Eq. (42) over x gives 



\/x' G X : P°°{x')T{x, x') = P°°{x')T{x', x) = P°° Yl T ( x '' x ) = P°°( x ')- (43) 
xex x<=x x<=x 

Hence, for reversible Markov chains, P°° is always a steady-state (stationary) distribution of ST . 
In the case where P°°(x) is uniform over X, i.e. P°°(x) = cste, the detailed balance condition 
reduces to 

Vx, x' G X : T(x, x) = T{x', x). (44) 

We thus find that if a Markov chain obeys detailed balance, there exists an invariant (stationary) 
distribution P°° over X. 



Appendix B A wrong way to define and use the Creutz algorithm 

A common error in using the Creutz algorithm is to define y l+1 as the first accepted move after 
y l with the condition Eq < $n{u 1+1 ) < Eq + AE, hereby discarding any unaccepted state in the 



expectation value of an observable A, c.f. Eq. ( 26 ) . We consider the Markov chain that corresponds 



to this wrong procedure and show that it does not verify detailed balance. 

Let us define more precisely the wrong algorithm. Given any configuration, we randomly pick 
a new configuration y l+l with probability T(x,x'). If Eq < <^n{u 1+1 ) < Eq + AE we accept the 
move. If the condition is not satisfied we keep picking at random new values for y l+1 until the 
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condition Eq < S^{x) < Eq + AE is fulfilled and then accept the move. This defines a Markov 
chain with transition probability 

R{x,x) = C(x')T(x,x'), (45) 

where C(x') = Y2x£Tm(E AE) T{x, x r ) is called the acceptance ratio (depending on x' only). 

The detailed balance condition would require \/x,x' £ Tn(E,AE) : R(x,x') = R{x',x). This 
would only hold when C(x') = C{x), but there is no reason for this to be true in general. 

We now recall how to properly empirically sample an observable A. The expectation value of 
observable A is computed through 

1 L 

<A(y)>= lim jY^Mv 1 )- (46) 

1=1 

In using the above defined wrong Creutz algorithm, the expectation value of A is calculated over 
Since the detailed balance condition does not hold for this Markov chain, one does not sample 
a stationary measure. 
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